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ABSTRACT 

We derive average flux corrections to the Model magnitudes of the Sloan Digital Sky 
Snrvey (SDSS) galaxies by stacking together mosaics of similar galaxies in bins of stel¬ 
lar mass and concentration. Extra flux is detected in the onter low snrface brightness 
part of the galaxies, leading to corrections ranging from 0.05 to 0.32 mag for the high¬ 
est stellar mass galaxies. We apply these corrections to the MPA-JHU (Max-Planck 
Institute for Astrophysics - John Hopkins University) stellar masses for a complete 
sample of half a million galaxies from the SDSS snrvey to derive a corrected galaxy 
stellar mass fnnction at z = 0.1 in the stellar mass range 9.5 < \og{M^,/M q) < 12.0. 
We And that the flux corrections and the use of the MPA-JHU stellar masses have a 
significant impact on the massive end of the stellar mass function, making the slope 
significantly shallower than that estimated by Li & White (2009), but steeper than 
derived by Bernardi et al. (2013). This corresponds to a mean comoving stellar mass 
density of galaxies with stellar masses log(M*/M 0 ) ^ 11.0 that is a factor of 3.36 
larger than the estimate by Li & White (2009), but is 43% smaller than reported by 
Bernardi et al. (2013). 
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1 INTRODUCTION 

The stellar mass function of galaxies is a basic probe of 
galaxy formation and evolution enabled by large redshift 
surveys. In recent years, major advances have been made by 
large redshift surveys, such as the 2dF Galaxy Redshift Sur¬ 
vey and the Sloan Digital Sky Survey (SDSS), in estimating 
the stellar mass function in the low-redshift Universe (Cole 
|et al.|2doT Bell et al.|2003 Blanton et al.|2003| . For exam¬ 
ple, Li & White 1 2009| have used a uniform sample of almost 


half a million galaxies from SDSS DR7 to derive the stellar 
mass function at z = 0.1. This has been complemented by 
the effort of the Galaxy and Mass Assembly Survey (GAMA 


Baldry et al. 20121, which has accurately constrained the 


faint end slope of the stellar mass function down to stellar 
masses ~ 10® Mq. 

The calculation of the stellar mass function hinges on 
the proper determination of the stellar mass of a galaxy, 
which in turn depends critically on the estimation of its to¬ 
tal flux in a given pass-band. Systematic differences in the 
estimation of the stellar mass of a galaxy may arise from 
different choices of the initial mass function (IMF) and the 
stellar mass-to-light ratio (M/L), as well as from different 
estimations of the galaxy total flux. Determining the flux 
accurately for a large number of galaxies in an all-sky survey 
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is a challenging task. In particular, quantifying the flux in 
the outer low surface brightness (LSB) regions of a galaxy 
has proven to be difficult and is still subject of much de¬ 
bate ( Bernard! et al.||2013 Simard et al.|[20TT l. These un¬ 
certainties mean that the slope at the massive end of the 
mass function is not very well determined. This has signifi¬ 
cant implications for several astrophysical problems, includ¬ 
ing halo occupation models, the mean baryon fraction in 
the Universe, X-ray and Sunavey-Zeldovich studies of high 
mass galaxies, and understanding the evolution of massive 
galaxies to high redshifts. 

Different approaches have been employed by SDSS in 
its photometric pipeline (PHOTO) to estimate the total flux 
of a galaxy. In addition to SDSS Petrosian magnitudes, 
two dimensional models (e.g. exponential or de Vaucouleurs) 
have been used to model the surface brightness distribution 
of galaxies (SDSS Model magnitudes). Further improvement 
has been provided by SDSS cModel magnitudes, for which 
fluxes are estimated as a linear combination of an exponen¬ 
tial and a de Vaucouleurs model. In recent years, several 
studies have tried to fit Sersic and multi-component models 


to the surface brightness distribution (Simard et al. 2011 


Lackner fc Gunn|2012| Bernardi et al.|2013l. 

Each of these approaches provides a progressively better 
estimate of the total flux of a galaxy, but they all suffer from 
the same intrinsic drawback, namely that the models are 
hts to the central, high signal-to-noise ratio (SNR) regions 
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of the galaxy and assumptions are required about the outer 
lower SNR (beyond Hr ~ 27 mag arcsec“^) part of the galaxy 
profile. Additionally, the total flux estimated through model 
htting can be biased in a number of ways. 

The biggest source of systematic bias in the flux deter¬ 
mination is related to the estimation of the sky background, 
especially for large nearby objects or those located in dense 
environments ( von der Linden 1|2007 Bernardi et al.|[2007 |. 
In principle, this can be overcome by considering extremely 
large fields of view. For example, considerable progress has 


been achieved by Blanton et al. (20111 by fitting the masked 


sky background for each SDSS scan with a smooth continu¬ 
ous function. 

However, even with improvements to the sky back¬ 
ground algorithm, one is still limited by the depth of the 
survey. The relatively short exposure time of SDSS (53.9 
secs) limits the accuracy of the background determination 
and subtraction. This in turn limits ones ability to distin¬ 
guish between the flux of the outer stellar halo and the sky 
background, leading to an over- or under-estimation of the 
total flux of a galaxy. In particular, multi-component model 
fitting of the main galaxy can lead to biased results. This 
may explain why recent attempts to trace the low SNR LSB 
part of a galaxy through fitting multi-component models to 
single SDSS photometric images have yielded divergent re¬ 
sults (Simard et al.||2011 Bernardi et al.|2013| Meert et al. 

^oTsl). 


Other sources of systematic error in determining the 
flux of a particular object are the procedures employed for 
deblending and masking, as well as the radial extent of the 
models used for the surface brightness fitting. Finally, in ad¬ 
dition to photometry, several other effects have a consider¬ 
able impact on the massive end of the stellar mass function, 
such as evolutionary corrections and fiber collisions (i.e. the 
fraction of galaxies not targeted for spectroscopy due to the 
fact that fibres cannot be positioned closer together than 55 
arcseconds on the SDSS plug plates). 

An alternative but viable approach to fitting models to 
individual images of galaxies, is to stack images of similar 
galaxies to quantify the average total amount of extra light 



galaxy stacks helps to reliably constrain the total amount 
of light especially in the LSB component. In addition, the 
background for stacked galaxies can be determined more 
accurately. This then provides a direct handle on the cor¬ 
rections to the Model magnitudes as a function of the stellar 
mass and galaxy type. 

In this paper, we attempt to derive flux corrections to 
the Model magnitudes and re-derive the galaxy stellar mass 
function at redshift z = 0.1 using MPA-JHU (Max-Planck 
Institute for Astrophysics & John Hopkins University) stel¬ 
lar masses (Kauffmann et aL]|2003 1 and the sample of Li & 


White (20091. We estimate corrections to the Model magni¬ 


tudes by stacking volume-limited samples in bins of stellar 
mass, concentration and model type. We also consider var¬ 
ious effects that may systematically bias the stellar mass 
function. 

In Section we define the samples used for deriving 
the corrections as well as the full sample used to derive the 


stellar mass function. In Section we derive the flux cor¬ 
rections to the Model magnitudes. In Sections and we 
derive the galaxy stellar mass function and the luminosity 
function respectively. In Sections and we summarise and 
discuss our results. Throughout this paper, we assume a flat 
ACDM cosmology, Dm = 0.25 and Qa ~ 0.75. We further 
assume a Hubble parameter h = 0.72 for the calculation of 
physical distance scales wherever necessary. 


2 SAMPLE SELECTION 


2.1 Sample for Calculating the Mass Function 


Following Li fc White] ( |2009 1, we select SDSS spectroscopic 
galaxies from the NYU-VAGC ( New York University - Value 
Added Catalogue) Q catalogue (Blanton et al. 2005[ ) with 
redshifts in the range 0.001 Si z ^ 0.5 and Petrosian r- 
band magnitudes in the range 12 Si rrirPet ^ 17.This 
gives us a total of 533442 galaxies, which are ideal for large 
scale structure studies. We further pruned the sample to 
523476 galaxies by retaining only those galaxies with a valid 
MPA-JHU stellar mass. We estimate the “effective” survey 
area to be 6570 deg^ (2.0084 steradians), by taking into ac¬ 
count the incompletness and the masked-out regions (due to 
bright stars) of the survey. 

For the stellar mass function, we use the stellar masses 
provided in the DR7 version of the MPA-JHU catalogu^ 
which assumes a universal Chabrier initial stellar mass func¬ 
tion ( Chabrier 2003). 

To derive the luminosity function, we use the r-band 
absolute Model magnitude (M,,o.i), corrected for evolution 
and K-corrected to its value at zo = 0.1 according to the 
following equation: 


M = m — DM{z) — K{z- zo) + QRz — zo ), 


( 1 ) 


where M is the absolute magnitude, DM{z) is the dis¬ 
tance modulus at redshift z, m the apparent magnitude, 
K{z\zo) is the A-corrections relative to a passband blue- 
shifted by Zo and the luminosity e-correction is parametrised 
linearly by Q^. The A-corrections were calculated using the 


code Kcorrect v4_3 (Blanton & Roweis 20071. In general. 


we assume a uniform luminosity evolutionary correction of 
Qr = 1.62 as derived by (Blanton et al.|[2003K 


2.2 Sample for Determining the Flnx Corrections 

To derive the corrections to the Model magnitudes, we stack 
volume-limited sub-samples of isolated galaxies defined from 
the parent sample in various ranges of stellar mass, concen¬ 
tration (A90/A50) and redshift (See Table [^. In each sub¬ 
sample, galaxies that were better fit by an exponential (Exp) 
or a de Vaucouleurs (deV) model by the SDSS pipeline (de¬ 
fined by comparing the likelihood values of the model fits 
from the SDSS PhotoObjAll database) were stacked sepa¬ 
rately. 

We select isolated galaxies by requiring that there are 
no brighter companions in the spectroscopic sample within 


^ Available at http://sdss.physics.nyu.edu/vagc/ . 

^ We also include the three survey strips in the Southern Cap. 
^ Available at http://www.mpa-garching.mpg.de/SDSS/DR7/ 
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Figure 1. Volume-limited samples used for stacking: Shaded con¬ 
tours show the distribution of the parent sample galaxies in the 
plane of stellar mass versus redshift. The seven coloured boxes 
indicate the redshift limits of the seven stellar mass sub-samples. 
These sub-samples are further divided by concentration (C) (not 
shown in the figure). The numbers in the coloured boxes indi¬ 
cate the fraction of centrals in these volume limited stellar mass 
sub-samples. 


We fit 2D axisymmetric models (single Sersic and don- 
ble Sersic models along with a constant background) using 
the Bayesian analysis described by D’Souza et al.| (2014l to 
individnal postage-stamp cntouts of the highest stellar mass 
and high-concentration galaxies (11.49 < log(M*/M 0 ) < 
11.69, 2.9 < C < 3.3 ) in the redshift range 0.14 < z < 0.18 
(covered by the sample G4 above - 38 galaxies) and 0.2 < 
2 < 0.4 (414 galaxies). The choice of the sample was moti¬ 
vated by the idea of testing the robustness of the Model mag¬ 
nitudes in the limits of high stellar mass and high redshift, 
where the relative contribution due to the sky background 
becomes increasingly significant. 

We compare our best fitting model with the Model mag¬ 
nitudes reported by the SDSS photo v5_4 pipeline. In Figure 
we plot a histogram of Mmodei — Mfu for each galaxy. The 
distribution is broad with a standard deviation of 0.25 mag¬ 
nitudes and is positively skewed. The median is shifted by 
is 0.03 magnitudes and the mean by 0.08 magnitudes. The 
large spread in the histogram arises from a degeneracy be¬ 
tween the best-fit model and the level of sky background. 
The results in Figure demonstrate that shallow single¬ 
exposure SDSS images are insufficient to accurately quantify 
the total amount of light in massive early-type galaxies to 
better than 0.25 mag. We also note that estimates of the 
total flux from a single SDSS image will also be affected by 
the deblending and masking algorithm]^ 

Because of the limitations in estimating total fluxes 
from single SDSS images, we have chosen to correct Model 
magnitudes using stacked images, where the increased 
signal-to-noise ratio better constrains both the model and 
the level of sky background. 


R ^ 1 Mpc (where R is the projected comoving separation) 
and \5z\ < 1000 kms“^. 

In Figurej^ the redshift limits of the various stellar mass 
sub-samples are shown projected along the plane of stellar 
mass versus redshift of the parent NYU-VAGC sample. The 
fraction of centrals in each stellar mass sub-sample are also 
shown. 


3 FLUX CORRECTIONS TO THE MODEL 
MAGNITUDES 


In this section, we derive corrections to the original SDSS 
Model magnitudes derived from the standard DR7 photo¬ 
metric pipeline (photo v5_4). We first demonstrate that the 
original SDSS Model magnitudes are biased, in agreement 
with other studies (Simard et al. |2011|[Bernardi et al.|201^ . 
We then proceed to derive corrections to the Model magni¬ 
tudes. 


3.2 Stacking images 


In order to derive the flux corrections to the Model mag¬ 
nitudes, we used the sky-subtracted SDSS Data Release 9 
images to create mosaics in the g, r and i bands centred 
on each galaxy in the sub-samples defined in Section |2.2[ 
The mosaics extend out to radii of 0.6 - 1 Mpc depending 
on the stellar mass and redshift range. We follow the stack¬ 
ing procedure outlined in D’Souza et al. (in preparation) 
and similar to that used by D’Souza et al. (20141. In short. 


each mosaic was deblended, masked, corrected for galactic 
extinction ( [Schlegel et al.|1998[ ), transformed to the highest 
redshift in that respective bin, rotated so that the major 
axis of each galaxy is aligned, and then stacked using the 
truncated-mean algorithm]^ The g- and the i- band images 
were only used to create the final mask along with the r-band 
images. Conservative masking was used. The final stacking 
was done using the masked and transformed r-band images. 


3.1 Systematic Biases in Model magnitudes 


The original Model are affected by two sources of systematic 


bias related to over-subtraction of the sky background (von 
der Linden 2007| Bernard! et al. 20071 and to simplistic 


choices of the model for the surface brightness profile of the 
galaxy (exponential or de Vaucouleurs). In this section, we 
allow for more complex models to describe the light profiles 
and we also allow the sky background to vary. 


3.3 Measuring the Total Flux of the Stacked 
Images 

Measuring the total integrated flux of a galaxy stack by 
fitting a model to its light distribution misses a fair amount 


In this paper, we follow the deblending and masking technique 
outlined in D’Souza et al. (in preparation). 

® For the truncated-mean stack, we removed 5 % of the extreme 
minimum and maximum values for each pixel. 
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Table 1. Volume-limited samples of isolated galaxies selected by stellar mass from the NYU-VAGC sample for the purpose of stacking 


Sample Stellar mass Concentration Redshift ^gal Exp ^gal deV 


Al 

9.69 < log(M./M 0 ) < 9.89 

1.7 

< 

C 

< 

2.5 

0.04 

< 

2 : 

< 

0.06 

797 

117 

A2 

9.69 < log(M./M 0 ) < 9.89 

2.5 

< 

C 

< 

3.3 

0.04 

< 

2 

< 

0.06 

66 

501 

B1 

9.89 < log(M*/M 0 ) < 10.09 

1.7 

< 

c 

< 

2.5 

0.05 

< 

2 : 

< 

0.07 

1028 

175 

B2 

9.89 

< log(Mt/M 0 ) 

< 

10.09 

1.7 

< 

c 

< 

2.5 

0.05 

< 

2 : 

< 

0.07 

83 

1111 

Cl 

10.09 

< 

log(M./A 40 ) 

< 

10.29 

1.7 

< 

c 

< 

2.1 

0.05 

< 

2 : 

< 

0.08 

638 

15 

C2 

10.09 

< 

log(M*/A 40 ) 

< 

10.29 

2.1 

< 

c 

< 

2.5 

0.05 

< 

2 : 

< 

0.08 

752 

308 

C3 

10.09 

< 

log(M*/M 0 ) 

< 

10.29 

2.5 

< 

c 

< 

2.9 

0.05 

< 

2 : 

< 

0.08 

121 

1499 

C4 

10.09 

< 

log(M./M 0 ) 

< 

10.29 

2.9 

< 

c 

< 

3.3 

0.05 

< 

2: 

< 

0.08 

2 

1071 

D1 

10.29 

< 

log(M*/A 40 ) 

< 

10.49 

1.7 

< 

c 

< 

2.1 

0.05 

< 

2: 

< 

0.09 

342 

38 

D2 

10.29 

< 

log(M*/M 0 ) 

< 

10.49 

2.1 

< 

c 

< 

2.5 

0.05 

< 

2: 

< 

0.09 

534 

535 

D3 

10.29 

< 

log(M./M 0 ) 

< 

10.49 

2.5 

< 

c 

< 

2.9 

0.05 

< 

2: 

< 

0.09 

89 

1468 

D4 

10.29 

< 

log(MfyA40) 

< 

10.49 

2.9 

< 

c 

< 

3.3 

0.05 

< 

2: 

< 

0.09 

1 

2153 

El 

10.49 

< 

log(M*/M 0 ) 

< 

10.69 

1.7 

< 

c 

< 

2.1 

0.06 

< 

2: 

< 

0.11 

239 

74 

E2 

10.49 

< 

log(M./M 0 ) 

< 

10.69 

2.1 

< 

c 

< 

2.5 

0.06 

< 

2: 

< 

0.11 

555 

1093 

E3 

10.49 

< 

log(MfyA40) 

< 

10.69 

2.5 

< 

c 

< 

2.9 

0.06 

< 

2: 

< 

0.11 

72 

1981 

E4 

10.49 

< 

log(M*/A 40 ) 

< 

10.69 

2.9 

< 

c 

< 

3.3 

0.06 

< 

2: 

< 

0.11 

- 

1867 

El 

10.69 

< 

log(M./M 0 ) 

< 

11.09 

1.7 

< 

c 

< 

2.1 

0.09 

< 

2: 

< 

0.13 

199 

264 

F 2 

10.69 

< 

log(M./M 0 ) 

< 

11.09 

2.1 

< 

c 

< 

2.5 

0.09 

< 

2: 

< 

0.13 

303 

1510 

F3 

10.69 

< 

log(M*/A 40 ) 

< 

11.09 

2.5 

< 

c 

< 

2.9 

0.09 

< 

2: 

< 

0.13 

76 

2919 

F4 

10.69 

< 

log(M*/A 40 ) 

< 

11.09 

2.9 

< 

c 

< 

3.3 

0.09 

< 

2: 

< 

0.13 

1 

4180 

Cl 

11.09 

< 

log(M./M 0 ) 

< 

11.69 

1.7 

< 

c 

< 

2.1 

0.14 

< 

2: 

< 

0.18 

6 

47 

C2 

11.09 

< 

log(M*/A 40 ) 

< 

11.69 

2.1 

< 

c 

< 

2.5 

0.14 

< 

2: 

< 

0.18 

11 

220 

C3 

11.09 

< 

log(M*/A 40 ) 

< 

11.69 

2.5 

< 

c 

< 

2.9 

0.14 

< 

2: 

< 

0.18 

15 

792 

C4 

11.09 

< 

log(M./M 0 ) 

< 

11.69 

2.9 

< 

c 

< 

3.3 

0.14 

< 

2: 

< 

0.18 

3 

2794 


of light due to the inability of the model to reproduce the 
bulge/disk component of the inner part of the galaxy. For 
example, a double Sersic model can miss upto 0.22 mag near 
the centre of a galaxy stack. On the other hand, “isophotal” 
magnitudes are unable to measure the LSB features of the 
galaxy stack, especially in the low S/N regime. 


In order to measure the total integrated light in each 
stack, we consider, therefore, a hybrid approach between 
a “model” and an “isophotal” magnitude. In particular, we 
first fit the large mosaics using two-dimensional axisymmet- 
ric double Sersic models with a flat background using the 
Bayesian analysis described by D’Souza et al. (20141. During 
this fit, the inner component is truncated outside a radius 
equal to 7 Re, while the outer component extends out to in¬ 
finity. Then, to the flux derived by the double Sersic model, 
we add the total residual flux (= Data —Model) within a cir¬ 
cular aperture of limiting radius Rum, defined as the radius 
at which the residual flux is maximised. The advantage of 
this hybird approach is that the double Sersic model mea¬ 
sures the slow decline of flux into the low S/N regime, while 
the sum of the residual and model fluxes reproduces the 
total flux in the high S/N part of a galaxy stack. 


In the next subsections, we explore the different factors 
that may bias our measurements of the total flux of the 
stack. 


3.3.1 Bias due to Inaccurate Sky Background Subtraction 


The sky residuals in the individual SDSS DR9 images are re¬ 
sponsible for some small amount of residual sky background 
in the final stacks (< 4 x 10“'* nanomaggies per pixel). We 
quantify this residual by adding a flat background compo¬ 
nent as a free parameter of our models (Section |3.3[ ). As show 
later in Section 3.3.2 this bias is minimal 0.01 mag). 

For the individual DR9 images, [Blanton et al.| ( |2011| ) 
quantified the spread in the sky background residuals to be 
a ~ 3.13 X 10“® nanomaggies per pixel. The median bias 
in the r-band magnitudes was estimated to be at the most 
0.1 mag independent of R50. In addition, higher stellar mass 
galaxies are found predominately at higher redshifts in the 
SDSS spectroscopic sample, limiting the bias in the flux of 
individual images caused by faulty sky-subtraction. 


3.3.2 Bias due to Models Used 

In order to test how the choice of model may affect the cor¬ 
rections to the total luminosity using the “hybrid” magni¬ 
tudes, we fit each of our galaxy stack using two-dimensional 
exponential, de Vaucouleurs, Sersic and double Sersic mod¬ 
els. Each model also includes a flat sky residual. By com¬ 
paring the evidences generated from the Bayesian fitting, 
we find that the double Sersic models are preferred by more 
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Figure 2. Bias in SDSS Model magnitudes: A histogram of the 
difference between the flux derived from our best fit models of 
high stellar mass high-concentration galaxies and the Model mag¬ 
nitudes from DR7 for galaxies with 11.49 < \og{Mt,/M q) < 
11.69, 2.9 < C < 3.3, 0.14 < 2 < 0.18 and 0.2 < 2 < 0.4. 
The green solid line indicates the fit to the skewed-normal distri¬ 
bution. 


than 10-a over the other models in all cases. The de Vau- 
couleur model gives the highest estimate of the total amount 
of light, followed by the single Sersic, the double Sersic and 
the exponential model respectively. Calculating the magni¬ 
tudes in the “hybrid” manner as described above gives very 
little difference in the total flux derived from different mod¬ 
els. 

Each model also yields different estimates of the residual 
background in the stacks. However, determining the back¬ 
ground level independently and keeping it fixed during the 
fitting process does not alter our estimates of the extra light 
(at the 0.01 mag level). This is due to the fact that the re¬ 
sults of the fitting are driven primarily by the inner high 
SNR part of the galaxy stack. 

We conclude that the combination of the depth of our 
stacked image and our “hybrid” magnitudes enables us to 
accurately constrain the total flux in the galaxy stack. Our 
outer models are not truncated, but instead extend out to 
inhnity. The difference between models which are truncated 
at 7 Re and models which instead extend out to infinity is 
at most 0.05 mag. 


3.4 Measuring the Flux Corrections 

For each stellar mass, concentration range and model fit type 
(exponential or de Vaucouleurs), we measure the average 
extra flux correction to the Model magnitudes as the differ¬ 
ence between the total integrated light in the stack and the 
median Model flux of the galaxies in the stack. The median 
Model magnitude was calculated by taking the median of the 


individual fluxes of galaxies in the stack. We find that the 
median Model magnitude is on average higher than the mean 
Model magnitudes. We use a two-dimensional-interpolation 
scheme to calculate the average extra light as a continuous 
function of stellar mass and concentration for each model 
type. These are shown in Figure]^ As can be seen, there is 
an extra light contribution from those galaxies which were 
fit by an exponential model both for high concentrations 
and for high stellar masses. The extra light correction from 
those galaxies fit by a de Vaucouleurs model comes predomi¬ 
nately from the massive, high concentration galaxies. On the 
other hand, the de Vaucouleurs model often over-estimates 
the flux of a galaxy for low concentration massive galaxies. 

We note that the large width of the stellar mass bins for 
the highest stellar mass galaxies may influence the correction 
derived in the stacking procedure. To account for this, we 
divide our highest mass sample, G4 (11.09 < log(M*/M 0 ) < 
11.69, 2.9 < C < 3.3, 0.14 < 2 < 0.18 ), into smaller mass 
bins of size 0.1 dex. We find that the relative corrections 
ranges from 0.23 to 0.31 mag, gradually increasing from the 
lowest to the highest stellar mass bin (see Figure]^. The 
mean correction derived by stacking the entire sample G4 is 
0.29 mag. 

For galaxies outside the the mass limits dehned in Sec- 
tion | 2 . 2 [ we extrapolate assuming the same mass corrections 
of the nearest defined mass bin. In particular, at the high 
mass end, there are 116 galaxies with stellar masses larger 
than log(M*/M 0 ) > 11.69, the highest stellar mass bin used 
above. For these galaxies, we assume the corrections to be 
the same as found for the highest stellar mass bin ( 0.31 
mag). 

Assuming a constant M/L for each galaxy, we calculate 
the extra mass for each galaxy in our main sample given its 
stellar mass, concentration and model type (by comparing 
the likelihoods of the Model fits from the SDSS database) 
as: 

log = -AMag/2.5 ( 2 ) 


4 THE STELLAR MASS FUNCTION OF 
GALAXIES 

4.1 Method 

We estimate the abundance of galaxies as a function of their 
stellar mass using the 1/Vmaa: method outlined by Li & 


White (20091. In combination with the depth and the large 
spectroscopic sample of SDSS, the 1/Vmax method provides 
an unbiased estimate of the stellar mass function and its nor¬ 
malisation.In Section [ 4 .2[ we demonstrate that the 1/Vmax 
estimator is unbiased against large scale structure at stellar 
masses of log(M*/M 0 ) ^ 9.5, which is the regime studied 
in this work. We limit ourselves to this regime since as esti¬ 


mated by Figure 4 of Baldry et al. (20081, all galaxies above 


stellar masses of log(M*/M 0 ) ^ 9.5, will be detected irre¬ 
spective of their central surface brightness. Moreover, our 
flux corrections begin from log(M*/M 0 ) 9.6 upwards. 

For each observed galaxy i, we define the quantity 
Zmax,i to be the maximum redshift at which the observed 
galaxy would satisfy the apparent magnitude limit of our 
sample nrirp^^ ^ 17.6. Evolutionary and K-corrections are 
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Figure 3. The flux corrections (AMag) as a function of stellar mass and concentration using an interpolation scheme for exponential 
(Exp) and de Vaucouleurs (DeV) fit galaxies. 





Similarly, we also define Zmin,i as the minimum redshift 
at which the galaxy would be present in our sample. Hence, 
Zmin,i is the maximum of the lower limit of the redshift slice 
and the solution to the equation: 

M, = - DM{Zmin) “ K{Zmin ) (4) 

This then allows us to calculate Vmax,i for the galaxy 
in question as the total co-moving volume of the survey be¬ 
tween Zmin,i and Zmax,i- The Stellar mass function can be 
then estimated as; 

= E(/" orm coll,i Vmax,i) ^ (5) 


where fnormcoH,i IS the normalised fiber collision factor de¬ 
fined below, and the sum extends over all sample galaxies 
with stellar masses in the range M*±0.5 5Mt. The error bars 
are estimated by taking into consideration both Poissonian 
and bootstrapping errors, as well as errors due to cosmic 
variance (See 4.21. 

We calculate the stellar mass function in the total red¬ 
shift range 0.001 ^ ^ 0.5 as well as in three redshift slices: 

0.001 s; 0 ^ 0.15, 0.15 ^ ^ 0.3 and 0.3 ^ ^ 0.5. 


Figure 4. The flux corrections (AMag) as a function of stellar 
mass for galaxies in the sample G4 (blue circles). We also indicate 
the uncertainty in the corrections by showing the flux corrections 
derived from the mean Model flux of the stacks (red circles). 


included when calculating Zmax,i- Hence, Zmax,i is the mini¬ 
mum of the upper limit of the redshift slice and the solution 
of the equation: 

Mi = - DM{Zmax) ^i^Zmax') Qei^Zmax Zi ) (3) 


4.2 Robustness of the IjVmax Estimator 

In this work, we estimate the abundance of galaxies using 
the XjVraax method. Given the large effective surface area 
(nearly 6570 degr^) and the depth of spectroscopic sample, 
the IjVraax method will be invariant to large-scale structure 
up to a limiting stellar mass. To test this, we divide our sam¬ 
ple into three independent but contiguous parts (Sample A, 
Sample B and Sample C split by right ascension), and calcu¬ 
late the standard deviation in the stellar mass function as a 
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Figure 5. Top panel: The stellar mass function for the complete 
sample as well as the three independent smaller samples A, B and 
C, using the MPA-JHU stellar masses for the full redshift range 
0.001 2 : s; 0.5. Bottom Panel: The standard deviation in the 

stellar mass function as a function of stellar mass. 


function of stellar mass. In the bottom panel of Flgure[^ the 
standard deviation is plotted as a function of stellar mass. 
The difference in the estimates of the stellar mass function 
due to the l/Vmax method from the three independent sam¬ 
ples Is less than 10% for stellar masses log^Mt/Mo) > 9.5. 
The standard deviation gives us also a handle on the er¬ 
rors in our estimates of the stellar mass function due to the 
cosmic variance. 


4.3 The Effect of Systematic and Random Errors 
on the SMF 

In calculating the stellar mass function, various systematic 
and random effects combine to affect the final result. We 
discuss each of these effects in turn in the following subsec¬ 
tions: 


4..3.1 MPA-JHU Stellar Masses and Extra light from 
Photometry 

The first source of systematic bias comes from the estimation 
of the stellar mass of individual galaxies. In this work, we use 
the MPA-JHU stellar masses to calculate the stellar mass 
function. This involves a change of flux (from Petrosian 


to Model magnitudes) and M/L ratio (from NYU-VAGC to 
MPA-JHU) relative to |Li &: White (20091. 

We find that the use of NYU-VAGC stellar masses 
based on the Model magnitudes rather than the Petrosian 
magnitudes introduces a shift beyond the knee of the stellar 
mass function towards a shallower slope at the higher mass 
end. This shift is then further increased when we switch to 
MPA-JHU stellar masses based on the Model magnitudes. 
The slope of the massive end of the mass function is shal¬ 
lower than that obtained by shifting the mass function de¬ 
rived from the NYU-VAGC stellar masses by 0.1 dex (See 
appendix of Li fc White|2009 AlogM, = 0.1). At a stellar 
mass of logM* ~ 11.5 Mq, this accounts for an increase in 
the stellar mass function by a total of 1.24 dex (a 0.57 dex 
Increase due to the change from Petrosian to Model mag¬ 
nitudes and a 0.67 dex increase due to the change from the 
NYU-VAGC to MPA-JHU M/L ratios). 

Our assumed M/L ratio affects our estimation of the 
stellar mass function. The use of the MPA-JHU stellar M/L 
ratios makes the slope at the massive end shallower than the 
NYU-VAGC M/L ratios. We note that the MPA-JHU M/L 
ratios are derived from models that include the possibility of 
complex star formation histories, whereas the NYU-VAGC 
assumes that red galaxies can be described by single stellar 
populations. Analysis of spectra of massive galaxies in the 
BOSS survey by Chen et al. (20131 indicates that the star 
formation histories of the most massive galaxies are charac¬ 
terised by episodic star formation histories. 

The extra flux derived from the photometry of stacked 
galaxies introduces a further shift, making the slope at the 
massive end of the stellar mass function even shallower. We 
find that this shift of the stellar mass function is indepen¬ 
dent of whether we apply the corrections only to the central 
galaxies, or to all the galaxies in the sample. Although a 
small difference is found at the knee of the mass function, 
both results are consistent with each other within the error 
bars. 


4.3.2 Fiber Collisions 

The second source of systematic bias is caused by fiber col¬ 
lisions. The NYU-VAGC catalogue lists the spectroscopic 
completeness fsp of each galaxy, defined as the fraction of 
photometrically defined target galaxies in the subarea for 
which usable spectra are obtained. The NYU-VAGC cata¬ 
logue calculates the average completeness for each of these 
subareas by taking into consideration overlapping plates. In 
the jargon of the NYU-VAGC catalogue, these subareas are 
called sectors, fsp contains information about the missing 
galaxies due to lack of fibers in dense regions, missing galax¬ 
ies due to spectroscopic failures, and missing galaxies due to 
fiber collisions. The average fsp for the sample defined above 
is 0.9146. However, fsp assumes that all galaxies with mea¬ 
sured spectra are randomly distributed within a sector, and 
hence cannot account for specific differences between high 
and low density regions in the same sector. In particular, due 
to fiber collisions, certain galaxies (e.g. satellite galaxies of 
large clusters found at high redshifts) will be preferentially 
missed. 

To account for fiber collisions, we define the fiber colli¬ 
sion fcou.i for each galaxy, as the fraction of photometrically 
defined target galaxies that fall within a area of 55” in ra- 
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dius. fcoU,i takes the values between 0.111 and 1.0 (that is, 
8 closest neighbours and no neighbours respectively). The 
average of fcoii,i over our whole sample is 0.93819. We nor¬ 
malise fcoii.i such that it’s average value is the same as that 
of fsp- fnorm coll,i = fuc* fcoii,i, where fac is defined as 
< fsp > / < fcoii > and takes the value 0.9749. The nor¬ 
malised fiber collision fnorm coii now has the general average 
properties of fsp, but can better account for fiber collisions. 

Weighting by normalised fiber collision maintains the 
normalisation of the stellar mass function at the low mass 
end and increases the mass function up to 22% at the high 
mass end. Li & White (20091 did not include fiber collisions 
in their derivation and we can only reproduce their stellar 
mass function by using Petrosian magnitudes and by ne¬ 
glecting the effect of fsp. 


4-3.3 Evolution Corrections 

The third main source of systematic error is related to the 
assumption about the passive evolution of galaxies both in 
their number density and luminosity. In order to construct 
a stellar mass function from a large redshift range ((0.001 
2 SS 0.5), we would need account for the passive evolution of 
galaxies using a so-called evolutionary correction. Assuming 
such a uniform evolutionary correction is problematic, since 
galaxy evolution is a function of galaxy type and cannot 
be described by a simple linear model. For example, star¬ 
forming galaxies will evolve more slowly in luminosity than 
early-type galaxies. 

In order to quantify the effects on the stellar mass 
function related to the assumptions about galaxy evolution, 
we consider two approaches. In the first approach, we as¬ 
sume a uniform evolutionary correction {Qr = 1.62), which 
would represent an upper limit for the evolution of early- 
type galaxies with high stellar masses and stellar popula¬ 
tions that evolve passively with time (i.e. in the absence of 
any mergers). In the second approach, we derive the stel¬ 
lar mass function without evolution in three redshift slices: 
0.001 ^ 2 Si 0.15, 0.15 ^ 2 Si 0.3 and 0.3 ^ 2 Si 0.5. 

In Figure we plot the stellar mass function derived 
using the MPA-JHU stellar masses, including a uniform evo¬ 
lutionary correction, accounting for fiber collisions and for 
the additional stellar mass corrections due to the extra light 
at large radii (red solid curve). In addition, we also indicate 
the mass function calculated in the three redshift slices men¬ 
tioned above, without evolution. As seen from Figure]^ the 
evolutionary correction has only a small effect on the stellar 
mass function (~ 10% at the massive end). This is related 
to the fact that the luminosity evolution is implicitly folded 
into the derivation of the M/L ratio. 


4-3.4 Uncertainty due to binning the data 

Another source of systematic bias is related to binning 
the data in calculating the mass function via the llVmax 
method. In particular, this introduces further uncertainty 
at the massive end of the mass function due to a combina¬ 
tion of the low number statistics and the steep slope of the 
mass function over this mass range. In order to quantify this 
uncertainty, we recalculate the mass function with different 
values for the bin sizes, from 0.05 dex to 0.4 dex. In particu¬ 
lar, larger bin sizes tends to bias the slope at the high mass 



Figure 6. The effect of evolutionary corrections on the stellar 
mass function: The red solid line shows the stellar mass func¬ 
tion calculated using MPA-JHU stellar masses corrected for miss¬ 
ing light from photometry of stacked galaxies, corrected for fiber 
collisions and with uniform evolutionary correction Q = 1.62 in 
the redshift range 0.001 Si z ^ 0.5. Also plotted are the stellar 
mass function without evolutionary corrections in redshift slices 
0.001 Si 2 Si 0.15, 0.15 Si 2 ^ 0.3 and 0.3 Si 2 si 0.5 dashed blue, 
green and violet lines. 


end of the mass function towards shallower values. Reduc¬ 
ing the bin size increases the steepness of the slope until a 
saturation limit of about 0.1 dex. The variation caused by 
changes in the bin size around the saturation limit is within 
the uncertainties derived by bootstrapping and within the 
Poissonian errors. Hence, we calculate the stellar mass func¬ 
tion in bins of 0.1 dex. 


4 . 3.5 Eddington Bias 

Another source of systematic bias in the stellar mass func¬ 
tion is caused by the random errors in the flux and M/L ra¬ 
tios of individual galaxies. Such an “Eddington” bias causes 
the stellar mass function to be higher in the low-number 
density part because of scattering from the lower stellar 
masses (higher number density). This becomes particularly 
acute because of the steepness of the stellar mass function 
at higher stellar masses. 

To correct for this bias, we assume a parametrized form 
for the stellar mass function. We convolve this function with 
a distribution of the uncertainties in the stellar mass. We 
then fit this convolved function to the binned values of 
the stellar mass function calculated from the data using a 
maximum-likelihood method. The best fit paramteric func¬ 
tion is thus our true stellar mass function corrected for the 
Eddington bias. For the parametric function, we assume a 
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Table 2. Parameters of a double Schechter function fit to the 
stellar mass function of SDSS galaxies. 


a 

(h^Mpc“^ logio 


logio M* 
(h-^Me) 


0.008579 -1.082 10.615 

0.000355 -1.120 10.995 


double Schechter function, given by 


'i'MdM 


'I'l 

M* 


g-M/M* 

-M/M* 



( 6 ) 


where dM is the number density of galaxies between M 
and M + AM. This provides a much better fit to the data 
relative to a single Schechter. We further assume that the 
uncertainties in the stellar mass are distributed normally in 
log,o(M*/M0). 

To estimate the uncertainties in the stellar mass, we 
first estimate the M/L uncertainties as a function of stel¬ 
lar mass from the MPA-JHU database. We find that the 
average uncertainty Alogjo(Af/I/) ranges from 0.08 to 0.1 
as a function of stellar mass. We then estimate the average 
uncertainty in the Model magnitude as a function of stel¬ 
lar mass. We find that the average uncertainty in the Model 
magnitude is ~ 0.02 mag across the stellar mass range con¬ 
sidered. Hence the M/L uncertainty is much larger than the 
flux uncertainty. 

We find that correcting the stellar mass function for the 
Eddington bias reduces it at the high mass end by as much 
as 0.48 dex. 


4.4 Results: Stellar Mass Function 


In Figurej^ we present our final estimate of the stellar mass 
function corrected for missing flux, fiber collisions, evolution 
and Eddington bias with that of the original Li & White 
(2009) in red and the Bernardi et al. (2013) (Sersic-Exp 
fits) stellar mass function in green. 

We provide a parametric representation of the stellar 
mass function for stellar masses greater than log(M,/M0) ^ 
9.5. The parameters of the double Schechter function are 
listed in Table An integration of our stellar mass func¬ 
tion for stellar masses greater than log(M*/M©) ^ 9.5 gives 
the mean comoving stellar mass density of the low redshift 
universe as </* = 3.7 ± 0.3 10®/iMpc“®. This amounts to 
a 35% increase in the mean comoving stellar mass density 
contributed from the same stellar mass range for the |Li fc| 
White (20091 stellar mass function. In particular, focussing 


on the high stellar mass end: the mean comoving stellar mass 
density of galaxies with stellar masses log(Mt/MQ) 11.0 
is a factor of 3.36 larger than the estimate by Li & White 
(2009), but is 43% smaller than reported by Bernardi et al. 
(2013). 



Figure 7. The stellar mass function: The blue solid line shows 
the stellar mass function calculated using MPA-JHU stellar 
masses, corrected for missing flux, fiber collisions, evolution and 
Eddington bias. The red dot-dashed line shows the original Li 
& White 2009 stellar mass function calculated using the NYU- 
VAGC stellar masses based on Petrosian magnitudes. Also shown 
are the Bernardi et al. (2013) stellar mass function values (green) 
based on Sersic-Exp fits to individual galaxies. The dashed ver¬ 
tical line indicates our lowest stellar mass limit above which the 
stellar mass function is not affected by surface brightness com- 
pletness issues. 


5 GALAXY LUMINOSITY FUNCTION 


Similar to the galaxy stellar mass function, we also calculate 
the galaxy luminosity function using the IjVmax method. 
However, more careful attention needs to be paid to the evo¬ 
lutionary corrections which affects the luminosity function 
not only via the derivation of Vmax, but also via the calcula¬ 
tion of a galaxy luminosity via equation We calculate the 
luminosity function using two approaches: in redshift slices 
(0.001 ^ 2 ^ 0.15, 0.15 ^ 2 ^ 0.3 and 0.3 sf 2 ^ 0.5) 
without evolution and using a uniform evolutionary correc¬ 
tion of Qr = 1.62. In Figure we present the results of 
Mo,I t- band luminosity function considering Model magni¬ 
tudes with photometric corrections from stacking, fiber col¬ 
lisions and evolutionary corrections in bins of 0.25 dex. We 
also indicate the luminosity funtion without evolution cor¬ 
rections in three redshift slices. A comparision of our results 


with those of Bernardi et al. (20131 would require a more 


careful treatment of luminosity evolution which is beyond 
the scope of this paper. 


6 SUMMARY 

In this paper, we have shown that stacking similar galax¬ 
ies together in volume-limited stellar mass and concentra¬ 
tion bins allows one to derive average flux corrections to the 
SDSS Model magnitudes. In particular, we find that these 
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Figure 8. The luminosity function; The Mo.i^ luminosity func¬ 
tion calculated with photometric corrections, fiber collisions and 
flux uncertainty in three redshift slices and assuming an uniform 
evolutionary correction of Qr = 1.62 

. We also show the corresponding luminosity function from 
Bernardi et al. (2013) Sersic-exponential fits. 


corrections range from 0.02 to 0.31 magnitude, depending 
on the stellar mass and concentration of the galaxy. 

We apply these corrections to the Model fluxes and 
re-derive the stellar mass function using MPA-JHU stellar 
masses, accounting for galaxy evolution corrections and fiber 
collisions. We find that the slope of the massive end of the 
stellar mass function is shallower than reported by liTlF 


White (20091, but much steeper than derived by Bernardi 
et al.| i20 T^ 


The biggest change in the slope at the massive end of 
the mass function comes from our adoption of the MPA-JHU 
stellar masses (as much as a 1.24 dex increase at logM, ~ 
11.5 Mq with respect to Li fc White|2009 |. This involves an 
Increase of 0.57 dex and 0.67 dex due to the changes in flux 
and M/L ratio respectively. The second major contributor 
is the bias caused by the uncertainty in M/L ratio and flux 
measurements of individual galaxies which accounts for a 
decrease of ~ 0.48 dex in the mass function at the massive 
end. Fiber collisions contributes to an increase of nearly 22% 
at the massive end. Galaxy evolution corrections accounts 
for a decrease of maximum 10% at the massive end of the 
mass function. 

We also derive the r-band galaxy luminosity function 
and obtain similar results. In particular, the biggest source 
of systematic uncertainty in the galaxy luminosity function 
is related to the model assumed for the galaxy evolution cor¬ 
rection. In this Paper, we use the evolution correction values 
derived by Blanton et al. (20031, which serves as an upper 
limit for galaxies at the bright end of the galaxy luminosity 
function. 


7 DISCUSSION 


The flux corrections to the SDSS Model magnitude and their 
respective uncertainties derived in this work by stacking mo¬ 
saics of similar galaxies in volume limited stellar mass and 
concentration bins are consistent with those presented by 


Simard et al. (2011). We find no evidence for the need of 


large flux corrections of the order of 0.5 magnitudes as pro¬ 


posed by Bernardi et al. (2013). 


Our results are also consistent with extremely deep 
imaging of nearby early-type galaxies, obtained with the 
MegaCam camera on the Canada-France-Hawaii Telescope 
which indicate that outer LSB light contributes 5 to 16 per¬ 
cent to a galaxy’s total luminosity ( Due et al.|2015 1. Stack¬ 
ing results for luminous red galaxies (average redshift of 
2; ~ 0.34) from |Tal &: van Dokkum| ( |2011[ ) also indicate 
that typical SDSS-depth Images miss about 20 percent of 
the total stellar light. 

A number of systematic differences could contribute to 
the discrepancy between our results and those by |Bernardi| 


et al. (2013). In the limit of low SNR, the determination of 


the sky background level can influence the measured flux of 
a galaxy derived from fitting models to the surface bright¬ 
ness distribution. The depth of an image limits ones ability 
to distinguish between the flux of the outer LSB features 
of the galaxy and the sky background, especially for large 
stellar mass galaxies at higher redshifts. The use of multi- 
component models aggravates this problem. 

The simultaneous estimation of the model parameters 
and the sky background level may be prone to system¬ 
atic bias, since these are often degenerate with each other. 


Bernardi et al. (2013) use the PyMorph algorithm (based on 


GALFIT), which estimates the galaxy flux based on model fit¬ 
ting along with a simultaneous estimation of the sky back¬ 


ground. Meert et al. (2013) and Meert et al. (2015) have 


already highlighted the effect of a bias in the sky subtrac¬ 
tion on the total flux of a galaxy. On the other hand, SDSS 
Photo pipeline estimates the Model magnitudes by first in¬ 
dependently estimating and subtracting the local sky back¬ 
ground. A similar procedure is followed by [Simard et al.| 
(20111. In this work, we use the background subtracted im¬ 


ages provided with SDSS DR9 to derive the flux corrections. 
In addition, the depth of our stacked images allows us to ac¬ 
curately determine the residual sky background. 

Estimating the total flux of a galaxy Is dependent on 
the exact procedure used for deblending and masking ( see 
[Blanton et al.pOllj and [Simard et al.||201ip . In particular, 
the amount of masking employed has a substantial effect on 
the amount of flux that is derived for a specific galaxy. In 
this Paper, we use the conservative masking described by 


D’Souza et al. (2015), which involves using multiple runs of 
SExtractor (Bertin fc Arnouts|1996 1. 


Guo et al. (2010) calculated the stellar mass function 


using the NYU-VAGC stellar M/L ratios and Model mag¬ 
nitudes using the methodology of Li & White (2009). The 


stellar mass function derived here has a large shift and shal¬ 


lower slope than Guo et al. (2010), owing primarily to the use 


of the MPA-JHU stellar masses and the flux corrections to 
the Model magnitudes. The results of our work will affect the 
majority of recent halo occupation and abundance matching 
studies (e.g. Moster et al. 2013) that use the measurements 
of the stellar mass function from Guo et al. (2010|. 
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Finally, we comment that the majority of studies of the 
evolution of the massive end of the stellar mass function 
have found suprisingly little change out to 2 ~ 1 (Maraston 
et al 2013; Moustakas et al 2013; Fritz et al 2014). The 
co-moving number density of galaxies with stellar masses 
greater than IO^^Mq has apparently remained constant over 
the past 9 Gyr, calling into question the late build-up of 
these systems through mergers and accretion. Our work has 
shown that a significant fraction of the mass of these systems 
may be “hiding” in low surface brightness outer components 
that are systematically missed by conventional photometric 
extraction software. Accurately quantifying the evolution of 
the stellar mass in these halos will be an important challenge 
for next generation deep imaging surveys. 
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